Rank abundance relations in evolutionary dynamics of random replicators 
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We present a non-equilibrium statistical mechanics description of rank abundance relations (RAR) 
in random community models of ecology. Specifically, we study a multi-species replicator system 
with quenched random interaction matrices. We here consider symmetric interactions as well as 
asymmetric and anti-symmetric cases. RARs are obtained analytically via a generating functional 
analysis, describing fixed-point states of the system in terms of a small set of order parameters, 
and in dependence on the symmetry or otherwise of interactions and on the productivity of the 
community. Our work is an extension of Tokita [Phys. Rev. Lett. 93 178102 (2004)], where the 
case of symmetric interactions was considered within an equilibrium setup. The species abundance 
distribution in our model come out as truncated normal distributions or transformations thereof 
and, in some case, are similar to left-skewed distributions observed in ecology. We also discuss the 
interaction structure of the resulting food- web of stable species at stationarity, cases of heterogeneous 
co-operation pressures as well as effects of finite system size and of higher-order interactions. 

PACS numbers: PACS 



I. INTRODUCTION 

Understanding the relationship between complexity 
and stability is a fundamental and controversial prob- 
lem in ecology [lj. Before the 1970s the proposition that 
highly complex communities are more stable than simple 
ones was widely supported @, However, this early 
intuitive idea was challenged by theorists in the 1970s, 
who discussed the stability of a community of species in- 
teracting randomly f4j. In particular, the applications of 
random matrix theory rigorously revealed that the sta- 
bility of a community strongly depends on complexity, 
e.g. diversity and statistical properties of the interaction 
matrix, such as variance and connectivity, and complex- 
ity tends to destabilize community dynamics Since 
then, many mathematical ecologists have studied random 
community models to explain the apparent contradiction 
between the complexity of real- world ecosystems and the 
results of these mathematical studies [H, 0] . Recent the- 
oretical developments, for example, have discovered sta- 
bilizing factors of random community models: compe- 
tition [8[ and antisymmetric prey-predator relationships 
Q. Empirical and theoretical works also suggested im- 
portance of omnivory (higher connectance) [IOL 1 1 11 ] and 
weak interactions [l2| for stability. 

If the relative abundances of the species in a commu- 
nity are measured, inevitably a small number of very 
common species will be identified (i.e. species with 
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a high abundance), along with some rare species and 
more numerous species of varying intermediate degrees 
of rareness. Clarifying the mechanisms underlying these 
rank abundance relations (RAR) (the relations between 
abundance and the number of species possessing that 
abundance) are clearly another fundamental problem of 
ecology HEEil. In conservation biology as well, knowl- 
edge of RAR helps one to predict the likelihood of pop- 
ulation persistence and community stability in face of 
global change. Various models have been applied to 
ecosystem communities [H, [H, 03, d, Oil and, in spe- 
cial, the recent progress of the theory of 'neutral' mod- 
els [13, [H [H H [H have aroused constructive 
discussions on theoretical predictions and the experimen- 
tal studies on RAR. As the neutral models mainly cover 
ecosystem communities where species compete for niches 
on a single trophic level like a tropical forest or a coral 
reaf, the models have left the more complex systems a 
mystery. Such systems occur on multiple trophic levels 
and include complex interactions, such as prey-predator 
relationships, mutualism, competition, and detritus food 
chains. Although RAR are observed universally in na- 
ture, their essential parameters have not been fully clar- 
ified. 

As a step to explore RAR theoretically, in this paper, 
typical rank abundance relations are derived using a ran- 
dom community model with few parameters such as the 
level of symmetry of interaction matrix and co-operation 
pressure or productivity. While random community mod- 
els can be criticized for a lack of immediate realism, they 
have the advantage of being exactly solvable by analyti- 
cal techniques. Random replicator systems have for ex- 
ample been considered as solvable models of interacting 
species in [27l |2{|. In particlar, species abundance 
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distributions of random replicator models with symmet- 
ric couplings have been computed in [29j using methods 
from equilibrium statistical mechanics. Such static ap- 
proaches are limited to cases of symmetric couplings be- 
tween species, in particular the presence of predator-prey 
pairs (for which interactions are highly asymmetric) can 
not be taken into account in such equilibrium approaches. 
In order to remedy these shortcomings, we here take a 
different dynamical approach, allowing for an extension 
to systems with an arbitrary proportion of predator-prey 
pairs. To this end we employ methods different from 
those of [29( and focus on an approach based on dynam- 
ical generating functionals and path integrals. 

It is interesting to note that stochastic models of com- 
plex dynamically assembled food- webs [3(| [HJ , which is 
from a simple dynamics governed by generalized birth 
and death events, derive reasonable species abundance 
distributions, in good agreement with real data. In such 
models the multi-species dynamics is effectively reduced 
to that of a representative species, subject to a 'mean 
field' interaction with the remaining system. In a simi- 
lar fashion our approach reduces the evolution of species 
randomly coupled via quenched interactions to a 'one- 
species' effective process as well (albeit a non-Markovian 
one). This mapping leads to an exact solution in the 
thermodynamic limit of infinite system size. For the 
stochastic approach, the model has the randomness with 
some probability. On the other hand, for generation func- 
tional, it gives the fixed randomness in the deterministic 
time evolution. Apart from providing a starting point 
for more realistic modifications of the present model, our 
analysis can hence, to a certain de gre e, be seen as com- 
plementary to the approach of [30l. l3l|. 

In the context of statistical mechanics another inter- 
esting point of the present model is that the replicator 
dynamics with asymmetric random interactions shows a 
non-equilibrium phase transition, i.e. two phases with 
qualitatively different behaviors are found (stable versus 
unstable). At the same time the replicator system does 
not exhibit a Lyapunov function, and is hence intrinsi- 
cally a non-equilibrium model without detailed balance. 
Further details can be found in the statistical mechanics 
literature [12, HH ■ In our system destabilization of a glob- 
ally fixed point solution and its bifurcation to limit cycle, 
heteroclinic cycle and potentially chaos is found when pa- 
rameters are varied. The random replicator model hence 
shows similarities, but also crucial differences compared 
with e.g. models of spin glasses [34], Hf| and neural net- 
work models [36l 13^ . It is hoped that the study of ran- 
dom replicator dynamics may hence contribute to the un- 
derstanding and classification of dynamical phase transi- 
tions in disordered systems. 

This paper is organized as follows: we will define the 
model in Sec. II and then discuss the statistical mechan- 
ics analysis based on a path-integral approach in Sec. III. 
In Sec. IV, we show results for pairwise interaction: a 
stability analysis, phase diagram, survival function, rank- 
abundance relations (RAR), the species abundance dis- 



tribution (SAD), finite size effects and structure of the 
resulting food web are discussed. We then turn to het- 
erogeneous co-operation pressure and higher-order inter- 
actions in subsequent sections V and VI, respectively. We 
summarize our results in Sec. VII. 



II. MODEL 

We here study the simplest system of random repli- 
cator subject to Gaussian interaction, and focus on 
the model originally proposed by Diederich and Opper 
[27l [28| . In conventional replicator dynamics, the sys- 
tem consists of N species, labeled by i = 1,...,N. 
The composition of the population of species at time 
t is then described by a concentration vector x(i) = 
(xi(t), . . . , Xiv(t)), where Xi(t) denotes the concentration 
of species i = l,...,N, and where Xi(t) = 1. The sys- 
tem evolves in time according to the following replicator 
equations [38[ 



Xi(t) 



= /,[x(t)] - u(t), 



(1) 



where /i[x] is the 'fitness' of species i at time t, and where 
v(t) denotes the mean fitness of species in the population. 
Hence species fitter than average increase in concentra- 
tion, whereas the weight of species less fit than average 
is reduced. 

We here take the fitnesses fi[x] to be frequency- 
dependent, i.e. they are functions of the vector x. Specif- 
ically we will assume, in the simplest setting, that 



/i[x] = -2uxi 



(2) 



i.e. that interaction between species is pairwise and char- 
acterized by the matrix elements Wj j . G eneralization to 
multi-species interaction is possible [22, Efj] , and will be 
discussed below. 

The matrix elements {iVij, Wji} (for any pair i < j) are 
chosen from a Gaussian ensemble. Specifically we choose 



= 0, 



; WijWji 



where 7TT denotes an average over the random couplings. 
w here characterizes the magnitude of the interaction, 
and r is a symmetry parameter and takes values T € 
[— 1 , 1] . For r = 1 the interaction between any pair of 
species i < j is fully symmetric, Wij — Wji. In this case 
no predator-prey pairs are found in the system. For T = 
and uiji are uncorrelated, the fraction of predator- 
prey pairs is hence 50 per cent. For T = — 1 all pairs 
of species are in predator-prey constellations, one here 



has 



Choosing intermediate values of T al- 



lows one to interpolate smoothly between these regimes. 
The ecologically most relevant setup corresponds to neg- 
ative values of T, describing prey-predator type interac- 
tion between species, rather than co-operation and direct 
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mutual competition. Diagonal terms in Eq. |2]) can be 
taken into account by writing wa = —2u, where u in the 
above setting denotes the so-called co-operation pressure 
[4l| . In an ecological context u takes mostly positive 
values. For u — > oo the ecosystem is found in a state 
of perfect co-operation and maximal diversity (with all 
species surviving and having equal concentrations). The 
essential parameter p = 2u can be termed as the pro- 
ductivity of a community in the sense of Lotka-Volterra 
equation (this will be explained in more detail in Sec. 
IV). Finally, in order to guarantee a well-defined ther- 
modynamic limit N — > oo, with which the statistical 
mechanics analysis of the model will be concerned, we 
re-scale the concentration vector by a factor of AT, and 
use the normalization iV _1 J2i x i(t) — !• Upon setting 
v(t) = iV _1 Xi(t) fi[x(t)] this normalization is con- 
served by the replicator dynamics ((T|). 

We will address the model by a combination of analyti- 
cal and computational methods. The statistical mechan- 
ics theory is described in the following sections, and its 
results will be compared against simulations in the subse- 
quent section. All simulations are here performed using 
the method described in [42j . This numerical scheme 
effectively amounts to a first-order forward integration 
with a dynamically adapted time-step. The latter is here 
necessary to avoid species concentrations to go negative 
in the discretized system. The dynamical time-stepping 
used in our simulations if typically of the order of 0.01 to 
0.1. 



III. STATISTICAL MECHANICS THEORY 



A. Path integral analysis 



are to be evaluated self-consistently as 

C(t,t') = (x(t)x(t')) ir , G( M ') = (|^) , (6) 

where (-)^ denotes an average over trajectories of the 
effective stochastic process The analysis then pro- 
ceeds by making a fixed-point ansatz, amounting to Q = 
C(t,t'). We also write X = J dtG(t) for the integrated 
response, and consider only ergodic states in which x re- 
mains finite. Restricting the analysis to asymptotically 
time-independent solutions of the effective process the 
following self-consistent equations for the resulting static 
order parameters Q, \ an d v (the fixed-point value of the 
average fitness) can then be derived similar to those re- 
ported in [H] 



M 

71 

QM 2 
A 

-M X 



[ Dz(A-z), 

J —OO 

/ Dz(A-z) 2 . 

J — OO 



Dz. 



(7) 
(8) 
(9) 



Here Dz — -A=e z ' 1 l 2 dz denotes the standard Gaus- 
sian measure, and one has A = w 2 Q, M = 2u + w 2 Tx 
and A = —vj\f\. We note that <f> = J_ Dz — 

h (l + erf (A/V2)) describes the fraction of surviving 
species. These equations are readily solved numerically, 
providing analytical predictions of the statistics of fixed- 
point solutions as functions of the model parameters w, T 
and u. 



The above system can be addressed by generating func- 
tional techniques originally devised in the theory of disor- 
dered systems [43[. It is also applied to linear evolution- 
ary dynamics model in [44j and can be adapted to the 
study of random Lotka-Volterra communities [HI]. We 
will not detail the mathematical steps here, as they have 
been reported in depth in the literature [28|, E3| ■ In the 
thermodynamic limit the system is found to be described 
by an effective single-species process of the form [28] 

x(t) = x(t)(-2ux(t)-T jf dt'G(t,t')x(t')-r)(t)-v(t)^ 

(4) 

(to denotes the time at which the dynamics is started). 
This process is non-Markovian in time, and subject to 
colored Gaussian noise rj(t), with temporal correlations 
given by 

( V (t) V (t')) = C(t,t'). (5) 

This colored noise is obtained from the interactions of 
randomness, which each trajectory has. C(t,t') and 
G(t, t') are the correlation and response functions, and 



B. Co-operation pressure and strength of 
interaction 

For reasons of completeness we re-iterate the phase 
behavior of the model as obtained by a linear stability 
analysis first reported in [28|. One here finds a stable 
region in which the fixed-point of the replicator dynam- 
ics is unique and locally attractive, separated from an 
unstable phase, as shown in Fig. Q] 

For r = — 1, the system is always found to be stable 
for any u > independently of w. At fixed w = 1, the 
onset of instability occurs at u c — v^2/4 and u c = V2/2 
for r = and T = 1 respectively. While the generat- 
ing functional approach is applicable for general symme- 
try parameter T, a static analysis based on the replica 
method is possible for symmetric couplings (r = 1). 
This has been carried out in [27], E^|. The replica ap- 
proach is here fundamentally different from ours, as it 
is only of a static (time-independent) nature. The er- 
godic fixed-point phase corresponds to a regime in the 
static analysis in which only one well-defined minimum 
of the Lyapunov function is found, corresponding to a 
so-called replica symmetric solution 35]. This solution 
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FIG. 1: (Color online) Phase diagram of the model with pair- 
wise interaction in the (w, u) plane for F — 1, 0, —1 from top 
to bottom. The system approaches a unique stable fixed point 
in the region above the respective lines, and remains unstable 
and non-ergodic below the phase boundaries. 



becomes unstable at the phase transition, referred to as 
a de Almeida-Thouless instability, coinciding with the 
location dynamical instability has been identified. For 
U c < u c (T = l,w = 1) = \/2/2 replica symmetry break- 
ing (RSB) occurs, i.e. the manifold of minima of the Lya- 
punov function becomes disconnected [35] . In conclusion, 
while the replica approach, requiring the existence of a 
Lyapunov function, is limited to the case r = 1, gen- 
erating functionals can be used to study the replicator 
system for any degree of asymmetry in the interaction 
matrix, as this approach requires only the knowledge of 
the dynamical equations of the system (the replicator 
equations), but is independent of the existence or oth- 
erwise of a quantity minimized by these dynamics. For 
the case of symmetric couplings, T ~ 1, the results from 
both methods coincide. 

Finally, since the behavior of the system can be seen 
to be qualitatively independent of the coupling strength 
w (which effectively re-scales the co-operation pressure), 
we will focus solely on w = 1 in the following. 



C. Species Abundance Distribution (SAD) 

Making a fixed point ansatz in the effective process (|3|) 
amounts to considering the time-independent solution of 
the effective species process of the form 28] 



x(z) 



M 



e 



(10) 



which represents the stochastic expression of the popu- 
lation in the stable state, z is here a static random vari- 
able drawn from a standard Gaussian distribution (%/Az 
reflects the single-particle noise r\{t) which becomes time- 
independent in the fixed-point regime). 0(-) is the step 
function. Note that, as mentioned above, only a frac- 
tion of species have positive concentrations at the fixed 



point, and that a complementary fraction of species dies 
out asymptotically. The distribution of concentrations x 
at the fixed point is thus a Gaussian cut-off at x = 
combined with a delta-peak at x = [42| . The so-called 
survival function 



a(x) 



lim 



1 

N 



"J2e(xi - x), 



(ii) 



denotes the fraction of species with a concentration 
strictly larger than x at the fixed point. The survival 
function, indicating the probability of a species having 
an abundance larger than x, is easily computed from (11C 
and is found as 



a{x) 




(12) 



in the thermodynamic limit. The fraction of survivors 
<f> as defined above is obtained as the special case <\> = 
a(x = 0). 

Using the cumulative distribution function C(x) 
1 — a(x) (denoting the probability for a species to have 
a concentration less than or equal to x) the abundance 
distribution for x > is given by 



F{x) 



dC(x) 
dx 



M 



: exp 



(A 



(13) 



A similar expression has been obtained for the case of 
symmetric couplings (r = 1) based on replica techniques 
in 29] . These earlier findings are found from our generat- 
ing functional analysis as a limiting case, so that genera- 
tion functional analysis contains the technique of replica 
method as mentioned above. 



IV. RESULTS FOR PAIRWISE INTERACTION 

A. Survival function 

We plot the survival functions a(x = 0) and a(x = 1) 
as a function of the co-operation pressure and for dif- 
ferent values of T in Fig. [5] As seen in the figure the 
diversity of the population (as measured for example by 
the number of surviving species) increases with larger co- 
operation pressure. The figure also demonstrates good 
agreement between numerical simulations and theoretical 
predictions for large values of the co-operation pressure 
u. In this phase the system is stable and ergodic and 
hence the fixed-point theory applies. Numerical simula- 
tions are performed using the discretization scheme de- 
scribed in [42} . Below a critical value u c (T) stability and 
crgodicity are lost (for T > —1), and the above theory 
can no longer be expected to be accurate, and system- 
atic deviations between theory and simulations may oc- 
cur. Still the qualitative agreement between theoretical 
lines, extended into the unstable phase, where they are 
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FIG. 2: (Color online) Survival functions a(0) and a(l) as 
functions of the co-operation pressure at fixed w — 1. Upper 
curves show a(0), lower curves a(l), with T = —1,0, 1 from 
top to bottom in each group. Lines are from theory (valid 
only above u c {T)), symbols from simulations of systems with 
N = 300 species, averaged over 20 samples. Surviving species 
in simulations are identified as species with Xi > 0.01 asymp- 
totically. 



technically no longer valid, is surprisingly good (RSB ef- 
fects have been seen to be weak in the low-u phase in 
previous studies). No unstable phase is present for fully 
anti-correlated couplings (r = —1) and non- negative co- 
operation pressure. 



B. Rank-abundance relations 

If the S = 4>N surviving species are re-labeled and or- 
dered according to their abundance in descending order, 
i.e. if X\ > X2 > ... > xs then a(x) can be understood 
as representing the species rank n according to 



a(x) 



N 



forx € [x n+1 ,x n ). 



(14) 



The function a(x) is a non- increasing monotonic func- 
tion, and can hence be inverted. The abundance x(n/N) 
of the n-th most abundant species can then be written 
as 



x{n/N) 



(n/N). 



(15) 



This representation is generally referred to as a 'rank 
abundance relations' (RAR) in the ecology literature. We 
find typical sigmoidal_patterns which have been observed 
in different regions [201 ] and with different species compo- 



sitions 47]. In general, for large value of u the RAR are 
broad and corresponds to RAR for a species-rich com- 
munity. Remarkably, the cross-over of the RAR patterns 
from low- to high- u is similar to the observed transition 
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FIG. 3: (Color online) Rank abundance relation for w = 1, 
r = — 1. Markers are from simulations. (N — 200, 20 sam- 
ples, 10000 iterations using the integration scheme of [42^1. 
lines from the fixed point theory. 
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FIG. 4: (Color online) Rank abundance relation for w = 1. 
F = 0. Markers are from simulations. (N = 200, 20 samples), 
lines from the fixed point theory, valid for u > u c — a/2/4, 
and of an approximate nature for u < u c . 



from low- to high productivity areas in real-world data, 
that is, comparing species-poor areas such as an alpine 
or polar region to a species-rich tropical rain forest [20j ] . 
The transition also corresponds to the secular variation 
of patterns observed in abandoned cultivated land [48| . 
This supports the contention that u is a maturity param- 
eter, as is suggested by an earlier evolutionary model in 
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FIG. 5: (Color online) Rank abundance relation for w = 1, 
r = 1. Markers are from simulations. (N = 200, 20 samples), 
lines from the fixed point theory, valid for u > u c — v2/2, 
and of an approximate nature for u < u c . 



FIG. 6: (Color online) F = — 1. The lines are from theory, 
u — 1.0, 0.8, 0.6, 0.4, 0.2 from top to bottom. The stable phase 
extends to all u > 0. Markers are from simulations. (N = 200, 
averages over 50 samples are taken). 



C. Species abundance distribution and Preston's 
octave plot 

Empirical data of species abundance have been taken 
for example in the studies of [5(| HH, H3, [Hj], and are 
normally presented as plots of 'species per octave'. I.e. 
species are grouped according to their abundance, and 
any species with abundance (number of individuals of 
that species present in the eco-system) in the interval 
of say [2", 2" +1 ) is subsumed in octave n (n being an 
integer). Log-normal distribution are then observed e.g. 
in [5 lL l52j . In order to depict the species abundance 
distributions in a manner similar to Preston's octave plot, 
we plot xF(x) versus x in a log scale following J2j|, see 
Figs. HE] and [Hill. 

Generally, we find that an increased co-operation pres- 
sure (equivalently an increased productivity, see below) 
larger u leads to 'octave plots' with small average and 
small variance. Species concentrations are here mostly 
found at a value of around x = 1 (in the limit of infinite 
co-operation pressure, u — ► oo, all species have equal con- 
centrations), and hence it is mostly the octave containing 
x = 1 which is populated. On the other hand, for smaller 
U, fewer species survive, and the variance in their concen- 
trations can be significant. This leads to octave plots of 
a large variance and a left-skewed form, similar to shapes 
observed e.g. in [iSll^l. In the fully asymmetric case 
r = — 1, see Fig. |5] all theoretical curves are in good 
agreement with results from simulations for all values of 
u. Here the theory is exact. In Figs. [7| and [5] however, 
corresponding to T = and T — 1 the theory is valid 
only for u > u c (T). Good agreement between analyt- 
ics and simulations is again observed. For u < u c the 
theory is at best of an approximative nature, and data 
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FIG. 7: (Color online) T = 0. The lines are from theory, u = 
1.0, 0.8, 0.6, 0.4, 0.2 from top to bottom. Stable phase contains 
u = 1.0,0.8,0.6,0.4. u = 0.2 is in the unstable phase where 
the theory applies only as an approximation (u c = v2/4), 
Markers are from simulations. (N = 200, 50 samples). 



from simulations appears much more prone to noise, and 
systematic deviations are observed from theoretical lines 
if they are continued into the unstable phase. Qualita- 
tively the theory is however able to capture the shape of 
the octave plots, in particular their left-skewness. 



D. Finite size effects 

Our theoretical analysis based on methods from statis- 
tical physics is mostly concerned with the limit of an infi- 
nite number of species in the ecosystem, N — > 00. This is 
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FIG. 8: (Color online) T = 1. The lines are from theory, u = 
1.0, 0.8, 0.6, 0.4 from top to bottom. Stable phase u = 1.0, 0.8 
and unstable phase u = 0.6,0.4 (u c = y/2/2). Markers are 
from simulations. (N = 200, 50 samples). 

of course for analytical convenience only, but can be ex- 
pected to be accurate also in the limit of large, but finite 
system size, as in real- world eco-networks. To study devi- 
ations from the exactly tractable infinite-size limit we dis- 
cuss simulation results of the species abundance distribu- 
tion of small systems in Fig. [5] One realizes that the dis- 
tribution becomes more left-skewed as the system size N 
is reduced, and that systematic deviations from the the- 
oretical lines emerge for systems smaller than about 100 
species. For smaller N, the amplitude of the peak gets 
larger. Note also that the largest possible concentration 
is limited by N (due to the normalization Xi = N), so 
that an effective upper cut-off is introduced for small sys- 
tems, and the distribution is skewed to the left. In nature 
it is impossible to obtain data for species with an infinite 
concentration, so that the part of the curve at small and 
intermediate concentrations seems most relevant. Sim- 
ulations indicate a trend toward more left-skewness at 
small system sizes. Unlike in other models of statistical 
physics at or near their phase transition points, we are 
unable to see fat-tailed broad species abundance distri- 
butions in the present model. 

E. Structure of the resulting food web 

TV-species replicator equations can in the context of 
ecology be shown to be equivalent to set of N — 1 coupled 
Lotka-Volterra (LV) equations [HI]. As discussed in [H, 
[54j the following transformation of variables 

yi = Xi/x M (t = l,2,...,JV) (16) 

n = W iM - W M M = W lM + P (17) 

bij = w i: j - w M j (18) 



FIG. 9: (Color online) F = 0, u = 0.4. The line is from theory, 
valid in the thermodynamic limit N — > 00. Parameters are 
chosen such that the system is in the stable phase, but close 
to the transition point of the infinite system. Markers are 
from simulations, N = 10, 20, 50, 200 respectively (averages 
over up to 10000 samples are taken for small system sizes). 

renders the replicator system studied in the previous sec- 
tions equivalent to Lotka-Volterra equations of the form 

= Vi \ Ti - Yl b vV^j ■ ( 19 ) 

The 'resource species' M £ {1, . . . , N} can here be cho- 
sen arbitrarily, note that one then has yu — 1 by con- 
struction, leading to an N — 1 dimensional system of LV 
equations. The ecological interspecies interactions bij are 
again of a Gaussian random form, but have different cor- 
relations than the couplings tUy of the original replicator 
system. For T = 1 and T = — 1 in particular, the b^ need 
not carry the symmetry (anti-symmetry respectively for 
r = —1) of the couplings wtj. The LV model describes 
an interaction network of species, where the interaction 
between any given pair of species (i ^ j) can be of 
a mutualistic type (bij and bji both positive) , of the com- 
petitive type (bij and bji both negative), or i and j can 
have a prey-predator relationship (one of the couplings 
positive, the other negative). These cases are summa- 
rized in Table HI The intraspecies interaction bu is given 
by bu = wu - w M i = -p - w M i = -r%. n is here the 
intrinsic growth rate of species i in the LV equations, 
and follows a Gaussian distribution of mean p and vari- 
ance 1/N. In particular, in finite systems, rj is positive 

with probability | ^1 + erf [yj N/2p)) . The parameter 

p(— 2u) can thus be interpreted as the 'productivity' of 
the community (the larger p the more species have posi- 
tive growth rate) . Note also that the average growth rate 

Y,i ~i is s iven b y P- 

In Figs. fTUl [TT1 and [T^l we depict the food webs in the 
stationary state of the replicator (or equivalently LV) dy- 
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Links 


green 


violet 


blue 


red 


(bij,bji) 


(+, +) 


(-,-) 


(+, ") 


(-. +) 


interaction 


mutual 


competitive 


i consumes j 


j consumes i 



TABLE I: Represented the interspecies interactions as colored 
links. Species are assumed to be ordered such that n > rj, 
for i < j. 



namics. Disks in these figures represent species, where 
species with a positive intrinsic growth rate (fj > 0) are 
shown as blue disks, and species with negative growth 
rate are depicted as red disks. Upon ordering surviv- 
ing species such that r\ > r 2 > ■ ■ ■ > > • • • > rs, 
the radius of the disk representing species i is chosen to 
be proportional to |Zo<7(ri)|/|/o<?|ri||. Note that the vari- 
ance of interaction strengths scales as I/TV in our model, 
i.e. Wij ~ 0(1 /vN), so for small u we can expect that 
|r*i| < 1 with large probability for any i (we have checked 
that \ri\ < 1 for all i for the data shown in Figs. [TUl fill 
and !12p . Since |^og|rj|| is monotonically decreasing func- 
tion of | r^l in the interval < |r;| < 1, larger blue disks 
hence mean larger productivity (i.e. fast growing species 
if interactions bij are switched off), and large red rep- 
resent large anti-productivity (i.e. species with quickly 
decaying concentration in the absence of interactions in 
the LV system). Links between species are shown in 
the figures only if the effective interaction exceeds a cer- 
tain threshold (i.e if max(|6y|, \bji\) > 0.6 * b max , where 
bmax =max(&jj) Vi, j). The thickness of each link in the 
Figures is in proportion to max(|6jj|, \bji\). 

The different types of interactions (see Table A} are 
represented by different colors: green links denote mutu- 
alistic interactions, violet competitive interactions, blue 
lines denote cases where a more productive species i ex- 
ploits a less productive one (j > i, assuming species are 
ordered such that r\ > r 2 > ... > r$) and red the reverse 
case of exploitation. 

In Fig. [TU] we depict a resulting food web for the case 
of symmetric interactions (T = 1), no red links are ob- 
served in this already reported in previous work 
On the other hand, one can see red links in Fig. 
ITT1 and Fig. [T3J For V = 1 the inter-species relation- 
ships are hence almost all mutualistic, i.e. there are no 
prey-predator type interactions in the equivalent Lotka- 
Volterra system. On the other hand, for T = — 1 the re- 
lationships are almost all of the prey-predator type and 
mutualistic enhancing interactions are found only very 
rarely in the Lotka-Volterra system. The case of uncor- 
related couplings in the replicator dynamics, T = 0, is 
an intermediate state. Finally, while we show the net- 
work topology only for small values of the co-operation 
pressure u in the figures, we note that with larger u, the 
network becomes more dense and of a more homogeneous 
structure. 




• • • 



FIG. 10: (Color) Network of interspecies interactions for w = 
1,T = 1,N = 100, w = 0.4. 




FIG. 11: (Color) Network of interspecies interactions for w = 
i,r = 0,N = 100, u = 0.2. 



V. SYSTEM WITH HETEROGENEOUS 
CO-OPERATION PRESSURE 

Heterogeneity between species is in the present model 
represented by the random interactions Wy. A second 
layer of diversity can be introduced, by making the co- 
operation pressure u species dependent, i.e. to use 

/,[x] = -2UiXi + ^ w ij x ] ( 20 ) 
ij 
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FIG. 12: (Color) Network of interspecies interactions for w 
l.r = -1,N = 50, u = 0.2. 



as the fitness of species i, where now Ui carries an ex- 
plicit index i and may be different from species to species. 
This model has been introduced and studied with gen- 
erating functional techniques in [5S| . In this section we 
will briefly discuss how adding heterogeneity of this type 
effects the distribution of surviving species, and will show 
how it can give rise to non-Gaussian abundance distribu- 
tions and how these can be computed from the statistical 
mechanics theory. Specifically we will draw the {ui} from 
a flat distribution over an interval [uq — cr, uq + er] , so that 
wo controls the mean co-operation pressure, and a > 
is variability over the ensemble of species. The gener- 
ating functional analysis is straightforward, but leads to 
an ensemble of effective species processes, one for each 
co-operation pressure present in the population. A fixed- 
point ansatz then leads to coupled equations for the static 
order parameters Q, \, A, expressed as integrals over the 
distribution of co-operation pressures, as reported in [55| . 
For x > the distribution of concentration of surviving 
species is then found as 



F(x) 



1 

2^ 



■ M(u) 
nil exp 



(A 



w 2 Q. 



(21) 

This is a su- 



where M{u) — 2u + w 2 T\, A 
per/position of cut-off Gaussians, with varying mean and 
variances, and may hence for sufficient width a of the 
distribution of co-operation pressures be of non-Gaussian 
shape. This is indeed observed in Fig. [T21 where we de- 
pict Fix) in a linear-log scale for various degrees of het- 
erogeneity in the co-operation pressures. For small values 
of the width a, the resulting function distribution F(x) is 
relatively close to being Gaussian, but can develop slowly 



FIG. 13: (Color online) Linear-log plot of distribution F(x) 
of surviving species for a system with heterogeneous co- 
operation pressures drawn from a flat distribution over [1 — 
a, 1 + a] where a = 0.1, 0.5, 0.75 from top to bottom at the 
maximum. Symbols are from simulations (T — 0, w — 1, 
N — 300 species, averages over 100 samples), solid lines from 
the generating functional fixed-point theory (note that for 
reasons of clarity we plot F(x) not xF(x) in contrast with 
other figures of previous sections). 
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FIG. 14: (Color online) Abundance distribution xF(x) for the 
system with 2-species and 3-species interaction. F = 0, the 
co-operation pressure is set as in indicated in the legend. 



decaying tails, and non-trivial kurtosis if the co-operation 
pressures become sufficiently variable across species. 
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VI. HIGHER-ORDER INTERACTION 

Up to now we have only considered the case of pairwise 
interaction between species. Generalization to higher- 
order interactions is possible and has been considered 
for example in [39l |40| . A random community model 
with p-body interaction between species can be defined 
as follows 



dt 



x l (t) = -Xi(t) 



2uXi(t) 



+ ^ , Ji 2 ,i3,...,i p X i2 {t) x i3 {fy ' ' ' X i v W 



(i 2 ,...,i P )eM. 



(p) 



-U(t) 



(22) 



with p a fixed integer and where = {(«2, • ■ ■ ,ip) ■ 

1 < 12 < 13 < • • • < i p < N\ %2, ■ ■ ■ , i p =f= i}. The cou- 
pling tensor is again assumed to be taken from Gaussian 
distribution with moments 



(Ji2,...,i p ) 2 



pi 



2NP- 1 ' 



T 2^- (23) 
We will consider p = 3 in the following. A generating 
functional and fixed-point analysis then leads to self- 
consistent equations 



Dz(A-z), (24) 
Dz(A-z) 2 , (25) 
Dz., (26) 



where Dz = —k=e z l 2 dz again denotes the standard 

V Z7T 



M 




7t 


r 

J — oo 


QM 2 


r 


A 


J — oo 






-M X = 


r 

J — oo 



Gaussian measure. These equations are very similar to 
the ones derived above for the case p = 2, differences 
are only to be found in the detailed expressions for the 
quantities M and A, which now read A = 2~i ^ 
'Id + 3TQx- We have A = —vj\f\ as before. 

Results for species abundance distribution of a replica- 
tor system with 3-species interaction are depicted in Fig. 
[H](for uncorrelated couplings, T = 0), and compared to 
the case p = 2 at otherwise unchanged parameters. For 
reasons of clarity we do not show results from numeri- 
cal simulations, even though we have performed numeri- 
cal tests in the ergodic stable phase and find reasonable 
agreement with the theoretical predictions. All other pa- 
rameters kept equal, a 3-body interaction appears to shift 
the peak of the distribution to the right, and to reduce 
its height, while increasing its width and left-skewness. 
Our findings thus suggest that higher-order interactions 
may add to the diversity of the ecological community, i.e. 
increase the variance of species concentrations at station- 
arity. 



VII. SUMMARY AND CONCLUDING 
REMARKS 



In summary we have presented a detailed discussion 
of species abundance relations resulting from the evolu- 
tionary dynamics of random replicator systems. Based 
on dynamical techniques of statistical mechanics of dis- 
ordered systems we have extended the work of [29|, [HJ 
to the case of asymmetric and anti-symmetric coupling 
matrices, and have also taken into account higher-order 
interaction modes and systems in which species are sub- 
ject to heterogeneous co-operation pressures. These sys- 
tems typically show a phase transition between a stable, 
ergodic regime and an unstable phase, in which the final 
state of the system depends on initial conditions. Based 
on a fixed-point ansatz the statistical mechanics theory 
is able to deliver exact analytical predictions for the re- 
sulting species abundance relations in the limit of infinite 
system-size, and computer simulations of the replicator 
dynamics are in perfect agreement with theoretical pre- 
dictions. The key findings of our analysis are the fol- 
lowing: (i) with larger co-operation pressure, regardless 
of inter-species interaction, the diversity of the popula- 
tion increases, (ii) we derive species-poor and species-rich 
RAR for symmetric interaction and species-rich RAR 
for asymmetric interaction, (iii) we find that the abun- 
dance distributions are typically similar to a lognormal 
distribution, and of a left-skewed type in our model, not 
too dissimilar from empirical data, (iv) visualizing the 
food-web structure of surviving species, and distinguish- 
ing between different types of pairwise species interac- 
tions gives insight into the stable relationship between 
species at stationarity, in particular symmetric interac- 
tions favor mutualistic relations, whereas anti-symmetric 
couplings tend to lead to one-sided exploitation of some 
species by others, (v) survival functions of systems with 
heterogeneous co-operation pressure can display highly 
non-Gaussian survival functions with long tails, (vi) in 
finite systems our theory is not applicable, and system- 
atic deviations are observed. In contrast with other dis- 
ordered systems SAD are not found to be fat-tailed or 
skewed to the right near the transition of the infinite-size 
model. 

The techniques we employ to study species abundance 
in random replicator systems are in the present con- 
text limited to fully connected random communities with 
Gaussian interactions. Extension to more realistic distri- 
butions of couplings may here be of interest, and simi- 
larly more realistic food-web topologies (see e.g. [56| or 
[HtJ and references therein) could be taken into account 
in future work. Methods from disordered systems the- 
ory can be adapted to those cases as well, and further 
studies would most likely be based on cavity methods 
or other tools used for finite-connectivity disordered sys- 
tems m, nil. 

There is currently also much interest in the relation- 
ship between deterministic models of population dynam- 
ics (defined through rate equations, e.g. the above repli- 
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cator dynamics) and stochastic individual-based models 
[oOl [6l| . It has here been seen that demographic stochas- 
ticity in models with a finite-number of individuals can 
induce behavior quite different from models based on rate 
equations It may hence be of interest to investigate finite 
microscopic individual-based analogues of random repli- 
cator systems (for example based on Moran dynamics) 
and to compare their dynamical behavior to that of the 
mean-field replicator system. Individual-based versions 
of systems with randomly drawn reaction rates have to 
our knowledge not been considered in the literature. This 
is indeed an interesting line of potential future work, al- 
though caution is appropriate when it comes to analytical 
approaches, as the randomness of interactions may make 
closed- form solutions of such models very difficult. 

It is hoped that our work may serve as a starting point 
for future studies in these directions, and that analysis of 
random community models of theoretical ecology based 
on methods from statistical mechanics may hence con- 
tribute to an understand of issues related to the diversity- 
stability debate as mentioned in the introduction. 
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Appendix 

The analysis of disordered systems by means of gen- 
erating functional is a useful and powerful method, es- 
pecially because it does not require the existence of a 
Lyapunov function, and is hence not limited to systems 
with symmetric interaction matrices. In this appendix we 
briefly outline the main mathematical steps and concepts 
of this technique. Further details can be found in a broad 
spectrum of sources in the literature [35j, |43j, l59j, |62| . 

The basic idea is to reduce a high-dimensional sys- 
tem with random couplings to an effective process for 



a respresentative (mean-field) particle. These processes 
are typically non-Markovian, even if the original system 
is Markovian, and subject to colored noise. If x(i) = 
(xi(i), . . . , xjv(i)) represents a trajectory of the micro- 
scopic system (subject to random interactions), then the 
starting point of the analysis is the dynamical partition 
function (or generating functional) 



exp 



(27) 



where (• • •) represents an average over all possible tra- 
jectories of the system. The dynamic partition function 
can hence be expressed as a path-integral over all such 
trajectories, and written in the form 



Z[ip] = J Dx <5(eq. of motion) e 



i£ t V(*K(*) 



(28) 



By 'equations of motion' we here mean the microscopic 
equations governing the dynamics, in our case the repli- 
cator equations Eq. (JTJ) , they contain the quenched dis- 
order of the problem (i.e. the random couplings). The 
analysis proceeds by writing the delta- functions in their 
Fourier representation by means of conjugate variables 
{xi(t)}, subsequently performing the average over the 
disorder, and then by introducing suitable macroscopic 
order parameters, such as e.g. the correlation function 
C(t,t') = N^ 1 J^i Xi(t)xi(t') and the response function 
G(t,t') = iN~ 1 J2i x i(t)%i(t)- I n the thermodynamic 
limit, N — > oo, an effective theory for C and G is then 
derived, expressed as a self-consistent problem involv- 
ing the above mentioned effective single-particle process 
in conjunction with self-consistent relations for correla- 
tion and response functions. As seen in Eqs. (|4l5p the 
effective process makes reference to C and G, and on 
the other hand these order parameters are to be com- 
puted self-consistently as averages over the ensemble of 
effective-particle trajectories $3§. 

For general systems the effective single-particle prob- 
lem can be addressed by suitable numerical schemes [63| . 
In the case of the replicator problem further analytical 
progress is possible based on the observation that the sys- 
tem attains a fixed point at sufficiently large co-operation 
pressure (28|. In this regime trajectories become effec- 
tively time-independent asymptotically, and further sim- 
plification is possible yielding Eqs. ([7][H]). Details of these 
steps can be found in [28[ and [40| . 
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